Hydrodynamics of confined colloidal fluids in two dimensions 
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We apply a hybrid Molecular Dynamics and mesoscopic simulation technique to study the dy- 
namics of two dimensional colloidal discs in confined geometries. We calculate the velocity auto- 
correlation functions, and observe the predicted t -1 long time hydrodynamic tail that characterizes 
unconfined fluids, as well as more complex oscillating behavior and negative tails for strongly con- 
fined geometries. Because the t -1 tail of the velocity autocorrelation function is cut off for longer 
times in finite systems, the related diffusion coefficient does not diverge, but instead depends loga- 
rithmically on the overall size of the system. 
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I. INTRODUCTION 

The role of hydrodynamics in two dimensions (2d) 
is considerably more complex than in three dimensions 
(3d). For example, when, in 1851, George Gabriel Stokes 
tried to extend his famous calculation of the low Reynolds 
(Re) number flow field around a sphere [l| to that of a 
cylinder he found that 

The pressure of the cylinder on the fluid con- 
tinually tends to increase the quantity of fluid 
which it carries with it, while the friction of 
the fluid at a distance from the cylinder con- 
tinually tends to diminish it. In the case of 
a sphere, these two causes eventually coun- 
teract each other, and the motion becomes 
uniform. But in the case of a cylinder, the 
increase in the quantity of fluid carried con- 
tinually gains on the decrease due to the fric- 
tion of the surrounding fluid, and the quan- 
tity carried increases indefinitely as the cylin- 
der moves on. 

so that there was no finite solution. This was later called 
the "Stokes Paradox". Experimental realizations of 2d 
systems are, of course, always embedded in one way or 
another in the 3d world. In a classic paper, Saffman Q 
demonstrated how taking into account the upper and 
lower boundaries on a 2d system solves the Stokes Para- 
dox because these boundaries open up a new channel for 
momentum flow out of the system. If the viscosity of the 
confining medium is 77', while the viscosity of the con- 
fined medium of height h is 77, then a new length scale 
emerges: 
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Stokes equations also cease to be valid at distances larger 
than L-^ e ~ v/U, where v is the kinematic viscosity 
and U the velocity of the fluid, because inertial forces 
must be taken into account. Although inertial terms 
also become relevant at similar length-scales in 3d, this 
fact doesn't need to be taken into account to obtain 
bounded solutions of the Stokes equations. For length 
scales L < min{Ls, Lj^ e } the total momentum in the 
2d layer is approximately conserved and Saffman showed 
that for a disk of radius R c and thickness h, the 2d dif- 
fusion coefficient for stick boundary conditions takes the 
following finite form Q: 
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beyond which the true 3d nature of the whole system 
needs to be taken into account. The zero Re number 



where fee is Boltzmann's constant, T is the temperature, 
and 7 = 0.5572 is Euler's constant. Note that in con- 
trast to the 3d form, where the diffusion coefficient only 
depends on fceT, R c and 77, here both the thickness of the 
film h and the viscosity of the boundary 77' enter into the 
expression for the diffusion coefficient. Eq. |T]) also im- 
plies that 2d hydrodynamics will be most evident when 
the confining boundary has a very low viscosity. 

Examples of experimental systems where 2d hydrody- 
namics are important include diffusion of protein and 
lipid molecules in biological membranes [HM IE 0, HI- 
Cicuta et al [§] recently directly measured the diffusion 
of liquid domains in giant unilamellar vesicles (GUVs) 
and found that the mean square displacement of the do- 
mains scaled logarithmically with their radius, in agree- 
ment with Saffmans prediction. 

Experiments on colloidal particles confined in a thin 
sheet of fluid (such as a soap film) have used video imag- 
ing [l(| and optical tweezers [ll| to explicitly demon- 
strate that the hydrodynamic interaction between the 
particles decays logarithmically with distance. These ef- 
fects can be understood from solving the 2d Stokes equa- 
tions and carefully taking into account the boundary con- 
ditions. Because the 3d boundary in these cases is air, 
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with a much smaller viscosity than the soap solution, 
Lg can be as large as 0.1m or more. The low Re num- 
bers typical of colloidal suspensions mean that L# e can 
be much larger than that, on the order of many meters. 
Furthermore, if the 2d systems under investigation has 
boundaries at a distance L <C ™in{£g, £j^ e } then the 
diffusion coefficient scales with system size as [l2j 

D ~ In [L/R c ] . (3) 

The goal of this paper is to use computer simulations 
to study the hydrodynamics of colloidal discs in confined 
geometries. We limit ourselves to two dimensions (2d), 
which has the advantage that simulations in 2d are faster 
than in 3d. The price we pay for this is that we must take 
into account some of the subtleties of 2d hydrodynamics 
described earlier, such as the finite size effects illustrated, 
for example, by Eq. ([3]). But these effects can also be ob- 
served in experiments on quasi-two dimensional systems, 
and are therefore interesting in their own right. 

We use a combination of Stochastic Rotation Dynam- 
ics (SRD) [HHIE! to describe the solvent, and Molec- 
ular Dynamics to solve the equations of motion for the 
colloids. Such a hybrid technique was first employed by 
Malevanets and Kapral [lal, a nd used to study colloidal 
sedimentation by ourselves [ItJ and by Hecht et al. [3 • 
We have recently completed an extensive study of this 
method to study the hydrodynamics of colloidal suspen- 
sions [15] . which we will call ref I, and we summarize 
some of the main points of the method in Section [TTJ 

Particle based methods like SRD (note that in the liter- 
ature this method is also sometimes called Multiple Par- 
ticle Collision Dynamics, see e.g. [H[). have the advan- 
tage that boundary conditions are very easy to imple- 
ment as external fields. This contrasts with traditional 
methods of computational fluid dynamics where bound- 
ary conditions are typically harder to implement. This 
suggest that methods like SRD may be ideally suited for 
the study of colloids in confined geometries. The rapid 
development of new methods to create microfludic sys- 
tems is also stimulating experimental studies on colloids 
in confined geometries [201 ] . For that reason, computer 
simulation techniques that can calculate the properties 
of colloids in narrow channels will become increasingly 
important. Another field of possible application includes 
flow in porous media [U, • 

We proceed as follows: In section [II] we describe the hy- 
brid Molecular Dynamics/SRD method we employ, and 
sketch out the key hydrodynamic parameters that govern 
the flow behavior. Section III describes simulations of a 
pure SRD fluid system in 2d, where we find that the ef- 
fects of hydrodynamic correlations are more pronounced 
than those found in 3d [19f . We also explore the impor- 
tant role of finite size effects. In Section IV we calcu- 
late the velocity autocorrelation function for colloids in 
2d, and show how confinement qualitatively affects their 
long-time behavior. In Section V, we analyze the dif- 
fusion coefficient for colloids in 2d, and show how the 
confinement effects seen for the velocity auto-correlation 



function are connected to the behavior of the diffusion co- 
efficient. We summarize our main conclusions in section 
VI. 



II. HYBRID MD-SRD COARSE-GRAINED 
SIMULATION METHOD 

To describe the hydrodynamic behavior of col- 
loids, induced by a background fluid of much smaller 
constituents, some form of coarse-graining is required. 
The hydrodynamics can be described by the Navier 
Stokes equations that coarse-grain the fluid within a 
continuum description. The downside of going directly 
through this route is that every time the colloids move, 
the boundary conditions on the differential equations 
change, making them computationally expensive to solve. 

An alternative to direct solution of the Navier Stokes 
equations is to use particle based techniques that exploit 
the fact that only a few conditions, such as (local) energy 
and momentum conservation, need to be satisfied to al- 
low the correct (thermo) hydrodynamics to emerge in the 
continuum limit. Simple particle collision rules, easily 
amenable to efficient computer simulation, can therefore 
be used. Boundary conditions (such as those imposed by 
colloids in suspension) are easily implemented as exter- 
nal fields. One of the first methods to exploit these ideas 
was direct simulation Monte Carlo (DSMC) method of 
Bird [H, (24|. The Lattice Boltzmann (LB) technique 
where a linearized and pre-averaged Boltzmann equation 
is discretized and solved on a lattice (25)], is a popular 
modern implementation of these ideas, and in particular 
has been extended by Ladd and others to model colloidal 
suspensions [H, M, H El, M, HI . 

In this paper we implement the SRD method first de- 
rived by Malevanets and Kapral [13j |. It resembles the 
Lowe- Anderson thermostat [32|, but has the advantage 
that transport coefficients have been analytically calcu- 
lated [II SI greatly facilitating its use. It is impor- 
tant to remember that for all these particle based meth- 
ods, the particles should not be viewed as some kind 
of composite supramolecular fluid units , but rather as 
coarse-grained Navier Stokes solvers (with noise in the 
case of SRD) [l|. 

An SRD fluid is modeled by N point particles of mass 
m, with positions and velocities Vj. The coarse grain- 
ing procedure consists of two steps, streaming and col- 
lision. During the streaming step, the positions of the 
fluid particles are updated via 

n(t + 5t c ) =Ti(t)+Vi(t)5t c . (4) 

In the collision step, the particles are split up into cells 
with sides of length a , and their velocities are rotated 
around an angle a with respect to the cell center of mass 
velocity, 

Vi(t + 8t c ) = v c . m ,i(t) + TZ t (a) [vi(t) - v c . Wl <(t)] (5) 
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FIG. 1: Top : Deviation of the simulated diffusion coefficient D s , from the random collision approximation D predicted by 
Eq. (|2f [) . as a function of the heavy particle mass. We simulated fluid particles in 2d for a square geometry with walls separated 
by a distance L = 32ao . Bottom : Deviation of the simulated diffusion coefficient D s , from the random collision approximation 
D , as a function of the particle mean free path A. Simulations were performed for values of A = 0.2, 0.4, 0.625, 0.8, f , f .2. 



where v c .m,i — Y] Z j (raVj)/y^ -m is the center of mass 
velocity of the particles the cell to which i belongs, Hi(a) 
is the cell rotational matrix and St c is the interval be- 
tween collisions. The purpose of this collision step is 
to transfer momentum between the fluid particles while 
conserving the energy and momentum of each cell. 

The fluid particles only interact with one another 
through the collision procedure. Direct interactions be- 
tween the solvent particles are not taken into account, so 
that the algorithm scales as 0(Af) with particle number. 
This is the main cause of the efficiency of simulations us- 
ing SRD. The carefully constructed rotation procedure 
can be be viewed as a coarse-graining of particle colli- 
sions over space and time. Mass, energy and momentum 
are conserved locally, so that on large enough length- 
scales the correct Navier Stokes hydrodynamics e merg es, 
as was shown explicitly by Malevanets & Kapral [13| . 

An advantage of SRD is that it can easily be coupled 
to a solute as first shown by Malevanets and Kapral [16| , 
and studied in detail in a recent paper by two of the 
present authors [lj| (ref I). If we wish to simulate the 
behavior of spherical colloids of mass M, they can be 
embbeded in a solvent using a Molecular Dynamics tech- 
nique. For the colloid-colloid interaction we use a stan- 



dard steeply repulsive potential of the form: 



Vcc{r) = 




while the interaction between the colloid and the solvent 
is described by a similar, but less steep, potential: 




where a cc and a cs are the colloid-colloid and colloid- 
solvent collision diameters. We propagate the ens uing 
equations of motion with a Velocity Verlet algorithm [35| 
using a molecular dynamic time step At 

Rt(t + M) = /WtJ + KMAf + ^At 2 (6) 

V,(W)=V,(,) t W. (7) 

where Rt and Vi are the position and velocity of the 
colloid, and F t the total force exerted on the colloid. 
Coupling the colloids in this way leads to slip bound- 
ary conditions. Stick boundary conditions can also be 
implemented [36| . but for qualitative behavior, we don't 
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expect there to be important differences. In parallel the 
velocities and positions of the SRD particles are streamed 
in the external potential given by the colloids and the ex- 
ternal walls and updated with the SRD rotation-collision 
step every time-step 5t c . 

To prevent spurious depletion forces, we set the inter- 
action range cr c / slightly below half the colloid diameter 
<7 cc /2 and include a small compensating potential for very 
short distances (when (3ip cc (r) > 2.5). For further details 
of how this procedure reproduces the correct equilibrium 
behavior see Ref I [lj| . 

The larger the ratio a cc /ao, the more accurately the 
hydrodynamic flow fields will be reproduced. Here we 
use cr cc /ao = 4.3, and a c f = 2ao, which was shown in ref 
I to reproduce the flow fields with small relative errors for 
a single sphere in a 3d flow. Other parameters choices 
taken from Ref I include e cc = e c j = 2.5fcsT for the 
colloids, and 7 — 5, a = ^ir for the SRD particle number 
density and rotation angle respectively. The time-steps 
for the MD and SRD step are set by slightly different 
physics (HI, and we chose At — 0.025io and St c — O.lio, 
where to = ao^J is the unit of time in our simulations. 

Coarse-graining methods like SRD are useful when 
they make the calculation of certain desired physical 
properties more efficient. To achieve this, compromises 
must be made (there is no such thing as a free lunch). 
For colloidal suspensions, for example, the Re number 
is typically very low, on the order of 10 -5 or less, and 
similarly the Mach number Ma = U/c Sl where U is a 
typical system velocity and c s is the velocity of sound, 
can be as small as 10 -10 . To achieve this in a parti- 
cle based simulation is extremely expensive. Resolving 
sound waves would mean that, since they travel much 
faster than colloidal particles, extremely small time-steps 
would be necessary in the simulation. Luckily even for 
Ma numbers as high 0.1 the hydrodynamics can be accu- 
rately approximated by incompressible hydrodynamics, 
so that one doesn't need to fulfill the physical condi- 
tion to obtain essentially the same physics. Similarly, 
for many applications, as long as the Re number is sig- 
nificantly lower than 1, the system can still be accurately 
described by the Stokes equations. A more detailed dis- 
cussion of these length-scales and hydrodynamic numbers 
can be found in ref I, and we will implicitly be making 
use of these arguments for the current work. 

A similar set of arguments can be made for the time- 
scales of a real colloidal fluid, compared to those found 
in our coarse-grained description. For example the kine- 
matic time, defined as t v = <j^. s /v, i.e. the time it 
takes a the vorticity to diffuse one colloidal radius, is 
of order 10 -6 s for a buoyant colloid of radius lfim sus- 
pended in water. For the same system, the diffusion time 
td = Vcc/D ~ 5s. Resolving these time-scales in one 
simulation would be very inefficient. In ref I we claim 
that successful coarse-graining techniques must telescope 
down the hierarchy of time-scales to a more manageable 
separations that are efficient for computational purposes. 



We argue that what is needed is not an exact represen- 
tation of all the time-scales of the physical system, but 
rather clear time-scale separation. For example, having 
t„ be only one or two orders of magnitude smaller than 
td can still lead to an accurate description of the desired 
physics. However, interpreting the results means taking 
this telescoping down of time-scales into account, and to 
do this properly, one has keep careful track of the physics 
involved. Expressing results as much as possible in terms 
of dimensionless units can facilitate this process [l5[ . 



III. DYNAMICS OF SOLVENT PARTICLES 

Before investigating the behavior of colloids in suspen- 
sion, we study a simpler problem of an SRD fluid confined 
to two dimensions. Much of this section will follow on an 
earlier comprehensive study by Ripoll et al. [19( in 3d, 
but here we focus on 2d. 

We begin by deriving an expression for the velocity 
autocorrelation function of the SRD particles, following 
similar steps to those found in ref. [l!| for 3d. The n th 
collision step of the SRD method can be rewritten as 

Vi{n5t c ) = Vi((n-l)St c ) + (ni(a)-X) 

x [vi((n-l)*t c )-v c . mji ((n-l)«t c )] (8) 

where X is the unit matrix, and t = n5t c the discretized 
time, with n the number of collision steps, St c the colli- 
sion interval and ~v c .m,i the cell center of mass velocity. 
The rotation matrix is defined in two dimensions as 

„ , N ( cos a ± sin a \ 
Tti(a) = . 

v ' V =F sin a cos a J 

such that the rotational average over any vector A be- 
comes 

((K(a)-l)A) = -(l-cosa)A = -( a {A). (9) 

If we now assume density fluctuations in each cell to be 
small, we can write (v c . m ,i(n6t c )) ~ ^(2j-'"v ). By 

multiplying each side by (vj(0)) and further assuming 
the velocity of colliding particles to be uncorrelated, we 
arrive at 

(vc.ro,i((n-l)ttc)v < (0)) ~ _L( Vi ((n- l)<yt c )v,(0)>. 

rwy 

(10) 

where 7 is the average number of solvent particles per 
cell. Substituting (p~0|) and ([9]) into (HJ), and rearrang- 
ing, we obtain an expression for the correlation of a fluid 
particle 

(vi(ntft c )vi(0)) = (1 - C«C)<Vi((n - l)5t c ) Vi (0)) (11) 

where £™ = 1 — 1/7. This expression shows that we can 
write the correlation at a certain time step in terms of 
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FIG. 2: The top two plots show the temporal evolution of the self diffusion coefficient of a fluid particle in a large box of size 
256ao x 256ao with periodic boundary conditions. The right plot shows that rather than saturate, the diffusion coefficient 
grows as D ~ \nt, as expected from theory. The bottom plot shows the hydrodynamic corrections to the diffusion coefficient 
D compared to to the random collision approximation expression Do given by Eq (|21|l for different box sizes L. As expected, 
these corrections shows a logarithmic growth with L/a . 



the previous time step, such that the normalized velocity 
autocorrelation function (VACF) is, 



(vt(n^ e )v i (0)) 
(vf(0)) 



(12) 



where £ = 1 — CaC™ is the decorrelation factor. The 
VACF, for reasons that will become apparent later, is 
the quantity of interest here and has the form 



(■ViinStcWO)) 



(13) 



A similar analysis can be performed for the case of a 
single heavy tracer particle of mass ml embedded in a 
solvent The total mass in a collision box is then 

(M + mj) such that the center of mass correlation is 
written as 



i(j*ft c )vi(0)> 



rwy + M 



'MnStMO)}. (14) 



By substituting (fT4f into (fTTj) . the decorrelation factor 
for a heavy tracer particle is found to be 



C = 1 - Ca 



rwy 



rwy + to' 



(15) 



The self diffusion constant D of a particle i is related 
to its mean square displacement via the Einstein rela- 
tion S3: 



D=lhn i([ r ^)-r 2 (0)] 2 ). 



(16) 



The position of a particle can be written explicitly in 
terms of discrete time-steps 



r i (*)=r i (0)+<5tc5^v i (Wt c ) 



(17) 



k=0 



so that 



n — 1 n — 1 

(k(t) - r,(0)] 2 ) = Stl J2 E (MjSQMkSQ). (18) 

j=0 k=0 



We note that combining the equation above with Eq. (jT6|) 
leads to the discrete form of the standard Green-Kubo 
expression for the diffusion coefficient as an integral over 
the velocity autocorrelation function. Manipulating the 
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sums, we find [3£ 

n— 1 n— 1 



EE 

j=0 k=0 



Vi(jSt c )vi(kSt c )} 



(19) 



71 — 1 n— 2 Ti—1 

E< v '0'^))+ 2 E E (viWt c )vi(k6t c )) 

3=0 3=0 k=j+l 

n-1 



2n- 



2^i(vi(0)v,((n-j)^ c )). 



Substituting the expression for the VACF derived earlier 
(fT3| into (fT9|) . we can write the diffusion coefficient in 
terms of its decorrelation factor £, 



D = lim — — — St r 



1 1 ^ ;r 

2 + nE^' 



2m [l-C. 

(20) 

Substituting Eq[l5] into Eq[20l results in the following 
dimensionless expressions for the self diffusion constant 
of a fluid and heavy tracer particle respectively 



D 



o 
Do 



= A 



Dn 



Am 
m' 



1 / rwy 
1 — cos a \ mj — 1 

1 (l + 



1 — cos a 
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(21) 
(22) 



Dq denotes the unit of diffusion and is expressed as 
a o/to = cto -y/fcsT/m and A is the dimensionless mean 
free path. It is a measure of the average distance the 
fluid particles travel in between collisions and has the 
form II 511 



x = s ± 

ao V m 



k B T 



to 



(23) 



These expressions for D make a key approximation, 
namely that collisions are always random, and that the 
particle velocities are uncorrelated. This neglects any hy- 
drodynamic effects. These expressions are thus expected 
to become more accurate if the mean-free path A be- 
comes larger so that the random collision approximation 
is expected to be a better description. Ripoll et al (l9j 
showed that in 3d, for their simulation parameters, the 
expression (|2"T|) for self-diffusion of an SRD particle began 
to show significant deviations from measured values when 
the mean-free path was smaller than 0.6. Similarly, they 
found that for smaller mean-free paths A = 0.1, these ex- 
pressions could underestimate the diffusion coefficient of 
a tagged heavier particle of mass M by as much as 75% 
for M > 10m. 

In Fig. [1] we analyze the self-diffusion coefficient of a 
tagged SRD particle as a function of mass and of mean 
free path for a square geometry with plates L = 32ao 
SRD cell widths wide. Similarly to Ripoll et al. [l!| we 
find deviations due to hydrodynamics, but in 2d these 



are much more pronounced. For example, as the mass in- 
creases, the hydrodynamic corrections to Eq|2"Tl saturate 
at a deviation of over 200% for larger masses. Similarly, 
we observe larger deviations as a function of mean free 
path than found in 3d. 

In contrast to the 3d results, for which finite size effects 
are not very strong, we expect that in 2d the effect of 
box size will be much more pronounced. To illustrate 
this, we carried out simulations in a much larger square 
box of width L — 256ao box sizes, now with periodic 
boundary conditions. These are shown in the top two 
plots of Fig. [21 We observe that the temporal diffusion 
coefficient, defined as 



D(t) 



{v(t')v(0)}dt', 



(24) 



continues to grow with time in a manner consistent with 
the expected scaling D ~ ln[t], as illustrated in the sec- 
ond top plot in Fig. [5J We expect the diffusion coeffi- 
cient to eventually saturate for this finite box size. But 
for an infinite box, we expect that D(t) will continue 
to grow indefinitely, a manifestation of the Stokes Para- 
dox. Similarly, for a fixed box of size L 2 , we expect that 
D ~ L/ao jl2| |. as discussed in the introduction, and this 
scaling is indeed observed in the bottom panel in Fig. [2] 



IV. VELOCITY AUTOCORRELATION 
FUNCTIONS OF COLLOIDAL PARTICLES 

Having worked out some properties of diffusing SRD 
particles, we now turn to the properties of colloidal par- 
ticles embedded in a solvent. 

If memory effects are ignored in a simple Langevin 
equation description of a spherical colloid of mass M, 
then the velocity autocorrelation function (VACF) of a 
colloidal particle can be calculated to be [39| 

(v(t)v(0)) = ^eM-t/h), (25) 

where the time — M/£ indicates how quickly particles 
forget their initial velocity. Its integral is related to the 
diffusion coefficient through the Einstein relation: 



D = 



(v(t)v{0))dt = 



(26) 



The Einstein relation is of course valid for any physical 
description of the VACF. 

Langevin approaches have traditionally been used for 
colloidal systems when hydrodynamics could be ignored. 
However, it is well known that hydrodynamic effects can 
have an important qualitative effect on the VACF. In 
their pioneering work, Alder and Wainwright (40| used 
MD simulations to demonstrate that the VACF (C(t)) 
of a tagged particle exhibits an algebraic decay at long 
times of the form t~ d / 2 , instead of the exponential form 
predicted by the Langevin equation. They showed that 
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t t/t 

v 

FIG. 3: Scaling of the velocity autocorrelation function : when VACF is normalized and plotted in terms of the reduced time 
t/t v , all the data collapse to the same curve. The VACF was originally measured for varying solvent densities 7. The system 
size is L — 32ao, which implies that \a — L/2a ca — 8. 



this behavior was a consequence of momentum conserva- 
tion, and therefore quite general. For colloidal particles 
in 3d, the diffusion coefficient is dominated by the con- 
tributions from this long time tail [la ], and we expect 
the same to be true in 2d. The correlation function for a 
colloid with slip boundary conditions can be calculated 
from kinetic theory (ioj : 



(v(t)v(0)) 



1 



dp J (4tt(L> + u)t) d l 2 ' 



(27) 



where d is the number of dimensions, and p the solvent 
density. This calculation predicts a t~ l power for the 
tail in 2 dimensions. That this should cause problems 
for the definition of D is evident from EqfSU because it 
implies that D diverges logarithmically with time. Note 
that similar behavior was seen for pure SRD particles in 
Fig. [2j where we found the scaling D(t) ~ ln[t}. For the 
colloids, we expect that the tail in the VACF will form on 
the timescale t„ = <J 2 s /v it takes the kinematic viscosity 
v to diffuse over the particle radius. 

FigO shows simulations run for a square box with a 
width L = 32gq. Eqf27] predicts that the tail should 
scale as ((v + D)pt)~ 1 . We tested this further by varying 
the number density 7 and simultaneously changing the 
density of the colloids so that they remain buoyant. For 
SRD the kinematic viscosity v depends only very weakly 



on 7 [l5l |34| for large values of 7 and keeping in mind 
that from equipartition 



c(o) = mo) 2 } = 



k B T 
M 



(28) 



it is not hard to show that the long time tails should all 
scale onto the same curve if time is scaled with t/t v . We 
show this explicitly in Fig. [3] for a fixed system size. 

At times shorter than the kinematic time, there is a 
contribution to the overall diffusion that comes from the 
local random collisions between the colloid and the sol- 
vent particles. This is typically dominant on time scales 
less than the sonic time t cs = v s /a cc over which collective 
modes can be generated [l5j]. We can calculate it using 
standard Enskog kinetic theory, and the ensuing Enskog 
friction coefficient has the following form [42| 



c 2d 



3%/2 



-cr cs 77T 



3/2 



k R T- 



mM 
n + M 



1/2 



(29) 



in two dimensions. Thus for very short times, the decay 
of the VACF is characterized by the Enskog time ts — 
M/£| d and it follows that 



(v(t)v(0)} 



k B T 
M 



exp(-i/£ B ). 



(30) 



because the collisions are essentially random. 
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FIG. 4: Normalized colloid VACFs simulated for increasing 
pipe width with Xs = = 7,24, with L = 28,96a . Sim- 
ulations results were scaled with the Enskog time is tE — 
1.0888to calculated from (I29|) . The dashed line represents the 
short time decay from (|30p . the dashed-dotted line represent 
the decay from the long time tail calculated from (|27p whilst 
the dotted line and the solid line denote the Langevin decay 
with the friction £ for the two respective box sizes. 
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FIG. 5: Log-log plot of the VACF for of a colloidal particle in 
boxes of size L = 32, 48, 500a (Xs = 8, 12, 125). The dashed 
line is the fit from Eg 1271 and serves as a guide to the eye. 



As shown in Fig. [U for short times, on the order of 
the Enskog time ig, the autocorrelation function shows 
clear exponential decay, in good agreement with Eq|30l 
The simulations shown are for two box sizes, and for 
short times, the VACF are independent of system size, 
as expected from Enskog theory. 



At longer times, Fig. [4] clearly shows the beginning of 
the long-time tail. The theoretical line we plot is from 
Eq. [57J an d fits remarkably well to the data. However, 
we note that there are some small deviations with system 
size at these longer times, which will be explained below. 

We also note that a direct comparison with the 
Langevin equation shows that for short times the 
Langevin equation overestimates the VACF, and that for 
longer times in underestimates the VACF for colloids. 
A more in depth discussion of this point can be found 
in appendix B of Ref I. In addition, in two dimensions, 
the Langevin equation (|25[) would predict an exponential 
form with different for different box sizes, because the 
diffusion coefficient changes with box size. By contrast, 
our results show that for short times the VACF is inde- 
pendent of box size. Clearly the Langevin equation does 
a poor job in capturing details of the colloidal VACF. 

While the simulations above are for fixed boundaries, 
it is also interesting to see what happens to the VACF 
when the confinement is more pronounced. In confined 
geometries, the particle induced flow fields should feel the 
presence of the walls. Bocquet and Barrat [43| showed 
that a sink in the decay of the long time tails should occur 
after an observation time on the order of t w = -4; = xt^u- 
This time is characteristic of how long it takes for the 
kinematic viscosity to reach the wall, when L/2 is the 
average distance to the wall. We illustrate the effect of 
the wall on the VACF in Fig. [5] for three different box 
sizes. For the two narrower boxes, the VACF clearly be- 
gins to drop below the i -1 power law but for the largest 
box, of size L = 500ao, we don't observe any deviation 
within our error bars. The sink in the tail for the Xs = 8 
simulation run begins at an observation times less than 
10t„, whereas the kinematic wall time in this instance is 
tw ~ 64i„. That may be because of other wall effects 
that kick in earlier for such a narrow box, or it may be 
that the cutoff in the algebraic decay is gradual and com- 
mences sooner than predicted by Bocquet and Barrat. 

In an important study, Hagen et al (4l| used Lattice 
Boltzmann simulations to investigate the VACF of a col- 
loidal particle between rigid walls and found qualitative 
deviations from the standard long-time tails. In particu- 
lar, for a sphere in a narrow enough cylinder, they found 
negative tails for the VACF C x (t) parallel to the walls 
that exhibited an algebraic decay like C x (t) ~ t~ 3 / 2 . 
Similarly, for a two dimensional disc between two plates 
they found C x (t) ~ t -3 / 2 , and for a three dimensional 
sphere between two plates they found C x (t) ~ t~ 2 . These 
exponents depend on the confinement, rather than on 
the overall dimension of the system. They explained 
the emergence of this negative tail with a simple mode- 
coupling theory that takes into account the fact that the 
sound wave generated by the colloid becomes diffusive. 
They further noticed that for slip walls, the normal be- 
havior was recovered, suggesting the origin of the neg- 
ative tail lies in the existence of velocity gradients near 
the wall. 

We performed simulations of colloidal discs in a pipe 
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FIG. 6: Normalized velocity autocorrelation function of a colloidal particle confined between two walls for increasing values of 
Xs ~ L/2a cs - Simulations here were performed for pipe widths L = 8, 10, 12, 24ao. The different plots denote the components 
of the normalized VACF parallel C x (t) (top) and perpendicular C y (i) (bottom) to the walls. Note that all plots are scaled 
with the kinematic time t v — a 2 a /u, and the time-scale for momentum to diffuse the distance between the walls t w = L 2 /4u. 
The minimum of the negative tails observed for C x (t) scale on top of each other when time is scaled with the sonic time t/t cs 
(insert in upper right panel). Similarly, the oscillating tails for C y (t) show the same period when scaled with t/tw (bottom 
right panel). 



of length 512ao with periodic boundaries in the x direc- 
tion and with two stick boundary condition walls at a 
reduced distance Xs — L/2u cs = 2,2.5,3,6, apart in the 
y direction and show the results in Fig. [5] We find a 
negative tail for C x (t), the VACF parallel to the plates, 
We find that the amplitude of the negative tail grows 
with increasing confinement. Furthermore, when time 
is scaled with t/t cs , the different correlation functions all 
show a minimum at about t « 3i cs , suggesting that sound 
waves are indeed the dominant cause of the negative tail, 
as suggested in [4l|- For the larger confinement shown 
here, Xs — 6, the VACF does show a rapid decay, but 
there does not seem to be a negative tail. This suggest 
that the diffusive sound wave mechanisms are still play- 
ing a part in the smallest (x s = 8) simulations of Fig. 
and may explain why the VACF decays on a shorter time 
t/tw than predicted by Bocquet and Barrat (43[. 

In the bottom two panels of Fig. [5] we observe oscilla- 
tory behavior for C y (t), the VACF perpendicular to the 
plates. This can be explained as follows: when a particle 
moves in the y direction towards the wall, it sets up a mo- 
mentum flow which can reflect off the wall and come back 
a time later to push the particle in the opposite direction. 



This effect should become more pronounced for stronger 
confinement, as we observe. To check this mechanism, 
we note that the walls introduce another length-scale, 

T 2 

tw = which is the time it takes vorticity to diffuse 
to the walls. If this reflection mechanism is at play, we 
would expect the period of the oscillations to reflect this 
time-scale. In the bottom right panel of Fig. [51 we ob- 
serve that when C y (t) is scaled with the time t/tw the 
oscillation minima indeed fall on top of each other, at 
least for sufficiently strong confinement. 

As discussed by Hagen et al [4l[, the C x (t) should ex- 
hibit a negative tail that scales like t~ for sufficiently 
strong confinement. In the upper plot of Fig. [7] we in- 
deed observe that the exponent is greater than t -1 , and 
consistent with i~2 ; as expected, although our data is 
not clean enough to confirm the exact exponent. Sim- 
ilarly, the final decay of the component C y (t) appears 
closer to t^ 1 than to 

Clearly confinement has an important effect on the 
long-time behavior of the VACF, and there may be fur- 
ther subtle effects that we have not yet been uncovered. 
It would be interesting, for example, to see how the an- 
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FIG. 7: Log-log plot of the VACF of a confined colloidal particle for boxes of size L = 8, 10, 12ao ( X« = 
The dashed line and the straight line have slopes i _1 and t~ 3 ^ 2 respectively. Only the positive values of 



2, 2.5, 3 respectively). 
-C(t) are plotted here. 



gular correlation functions, studied in ref. [36| with SRD 
for 3d stick boundary colloids in the bulk phase, would 
behave under confinement. However, for the calculation 
of long-time tails, methods like Lattice Boltzmann tech- 
niques used by Hagen et al. [HI], where noise does not 
play a big role, may be simpler and faster to use. 



DIFFUSION COEFFICIENT OF COLLOIDAL 
PARTICLES UNDER CONFINEMENT 



The Einstein relation (j2Uj) directly relates the VACF 
and the diffusion coefficient. We found that for short 
times, the VACF was well described by an Enskog 
form ([30]) that was largely independent of the bound- 
aries, and that at longer times it exhibited a long time 
tail that was much more sensitive to the boundaries. For 
strong confinement, the tail could even be negative or 
oscillatory, but for weak confinement, it appears to scale 
asC(t)~t- 1 . 

For an unbounded 2d system, the diffusion coefficient 
does not converge, instead its behavior with time can be 
approximated as: 



D 2d {t) 



(v(t)v(0))dt 



o 

k B T 
M 
k B T 

(2d 
SB 



exp(—t/tE)dt 
o 

miL 

87T77 L J *" 



k B T 

%-Kpvt 



(31) 



where we have assumed that the Enskog and hydrody- 
namic contributions to the VACF can be separated (this 
is not quite true) and moreover that the hydrodynamic 
tail does not kick until a time scale on the order of the 
kinematic time t v . We also assume that D <C v. 



A. Simulations in the 'bulk' 

In Fig. [8j we present the temporal evolution of the 
self diffusion coefficient of a colloid for a large box. We 
approximated colloids in the bulk by using a box of size 
L 2 = 256ao x 256ao with periodic boundary conditions. 
The plot shows results for solvent densities 7 = 5, 10, 50. 
On the time-scales of the simulation, we observe behavior 
consistent with D ~ ln[i], as expected from the t~ x tail 
of the VACF. In practice this would mean that D would 
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FIG. 8: Time evolution of the time-dependent diffusion coefficient of a buoyant colloid in two dimensions for SRD particle 
densities 7 = 5, 10, 50. The bottom two graphs show the diffusion coefficient normalized by C x (0) = ksT/M, with time scaled 
by t/t v . We expect the graphs to scale on to each other for longer times where the long-time tails dominate, (see also Fig. [3jl . 
The bottom right graph has a logarithmic scale and the dashed line has the slope ^f- and serves as a guide to the eye. 



grow indefinitely with time, and be unbounded, which is 
a manifestation of the Stokes Paradox. 



B. Simulations in confinement 

Whereas the diffusion coefficient of a two dimensional 
disc in the bulk appears to grow in an unbounded fashion 
with time, the diffusion coefficient for a confined fluid is 
expected to saturate at a finite value [ljj |43| . We showed 
in Figs [3] - [7] that the VACF is affected by the presence 
of walls, and no longer shows the t -1 behavior at very 
long times that would lead to a logarithmic divergence. 
As a result of the wall interaction, the diffusion will no 
longer diverge, but will plateau at a value determined by 
the distance to the wall. 

We tested this simple argument by simulating colloids 
under two different levels of confinement. The top panel 
of Fig. shows the integral of the velocity autocorrela- 
tion function plotted for colloids diffusing between par- 
allel plates a distance L = 32ao and L = 64ao apart 
respectively. For the smaller system, the temporal diffu- 
sion coefficient reaches a plateau at shorter times than is 
found for the larger system. 

To make these arguments more quantitative, we make 



the following approximation to the diffusion coefficient: 



D 2d (L/a cs ) 



(v(t)v(0))dt 



k B T k B T t 

— KT- H filth 

87T77 1 Jt » 

k B T t w 
D E + m ■ 



87777 t v 

k B T L 
In 

Anr] er cs 



(32) 



which indicates that the diffusion of a particle in confine- 
ment should scale with the log of the ratio of its radius 
to the pipe width. 

We performed simulations to check the validity of this 
simple scaling argument. The results are shown in Fig.[9l 
and can be accurately fitted to Eq. (|32f . 

While Eq|32] works very well for the larger boxes, 
it overestimates the diffusion coefficient for the smaller 
boxes. This is because the more complex wall effects 
shown in Fig [6] come into play so that the VACF no 
longer shows simple t _1 behavior assumed in Eq[32] We 
also note that in the smallest systems studied, the Enskog 
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FIG. 9: Top : The effect of confinement on diffusion. The two curves represent the temporal evolution of the self diffusion 
coefficient for colloids in varying confinement. The dashed line denotes colloids diffusing between plates a distance L = 32ao 
apart (xs = 8) whilst in the case of the solid line the plate separation is 64ao (Xs = 16). For the smaller pipe, the kinematic 
wall time is tw ~ 64i„. The diffusion coefficient begins to plateau much earlier than that because the effect of the wall on the 
VACF kicks in earlier than that. Note that only the x-component of the diffusion is plotted here. Bottom : Simulations were 
performed for 2d pipes of respective widths L = 8, 12, 16, 20, 24, 28, 32, 64, 96ao. The straight line is the fit from (|32fl and the 
circles are the simulation points. 



contribution is more than half of the overall diffusion (al- 
though this is not the reason for the deviation from the 
simple ln[X/(7 cs ] scaling). 

It was shown by Bungay and Brenner [44| using stan- 
dard methods of low Re number hydrodynamics, that the 
diffusion coefficient of a sphere in a 3d pipe of radius R p 
drops rapidly with R p /R c . Here we see from direct simu- 
lations of the VACF and the diffusion coefficient that the 
same behavior can be seen in 2d. It also drops as a func- 
tion of decreasing pipe width. It would be interesting to 
see if similar hydrodynamic arguments to those used by 
Bungay and Brenner [44[ could be used to explain the 
more rapid decrease of the diffusion coefficient observed 
in 2d for stronger confinement. 



VI. CONCLUSION 



strong confinement. We also observed the logarithmic 
dependence of the diffusion coefficient on system size, as 
originally predicted by Saffman [3[ for the lateral diffu- 
sion of a cylinder in a film. 

Although the Saffman result describes the motion of a 
disk of thickness h, and our simulation deals with disks, 
we can still map our results onto a real physical system 
by equating the diffusion coefficient measured in our sim- 
ulations to that measured in experiment. 

Through this study we have shown that SRD can be 
fruitfully used to simulate colloids in two dimensions. 
This suggests that it could easily be adapted for the study 
of other problems, such_as_ protein and lipid molecules in 

, liquid domains in gi- 



zn as protei 
biological membranes 0, H,H, 0, 
ant unilamellar vesicles (GUVs) 



, or colloids in a liquid 



film [n| HH , or various examples from microfluidics [2(| • 



We have applied the SRD simulation method to the 
study of the dynamics of two dimensional disks in con- 
fined geometries. We calculated the VACF for colloids 
and observed the predicted t~ x behavior as well as the 
more complex oscillating behavior and negative tails in 



Acknowledgments 

J.S. thanks Schlumberger Cambridge Research and 
IMPACT FARADAY for an EPSRC CASE studentship 
which supported this work. A.A.L. thanks the Royal 



13 



Society (London) and J.T.P. thanks the Netherlands Or- support. We thank E. Bock and I. Pagonabarraga for 
ganization for Scientific Research (NWO) for financial helpful conversations. 



G. G. Stokes, Trans. Cambridge Philos. Soc. 9, 8 (1851). 
Reprinted in Mathematical and Physical Papers, 2nd ed., 
Vol. 3 (Johnson Reprint Corp., New York, 1966). 
Cited in H. Lamb, Hydrodynamics, (CUP, Cambridge, 
1932), p 614. 

P.G. Saffman and M. Delbruck, Proc. Nat. Acad. Sci. 
72, 311 (1975); P.G. Saffmann, J. Fluid Mcch. 73, 593 
(1976). 

R. A. Cone, Nature New Biol. 236, 39, (1972). 

M. Poo, R. A. Cone Exp. Eye Res. 17 503, (1973). 

M. Poo, R. A. Cone Nature 247 438, (1974). 

M. Edidin, Annu. Rev. Biophys. Bioeng. 3, 179 (1974). 

R. Peters and R.J. Cherry, Proc. Natl. Acad. Sci. 79 4317 

(1982). 

P. Cicuta, S. L. Keller, S. L. Veatch; Journ. Phys. Chem. 
Lett B. Ill 3328 (2007). 

C. Cheung, Y.H. Hwant, X-l. Wu, and H.J. Choi Phys. 
Rev. Lett. 76, 2531 (1996). 

R. Di Leonardo, S. Keen, F. Ianni, J. Leach, M.J. Pad- 
gett, and G. Ruocco, Phys. Rev. E 78 031406 (2008). 
J. Happel and H. Brenner, Low Reynolds Number Hydro- 
dynamics (Springer, 1973). 
E. Guazzelli, Phys. Fluids 7, 12 (1995). 
2574 (1997). 

A. Malevanets and R. Kapral, J. Chem. Phys. 110, 8605 

(1999) . 

T. Ihle and D. M. Kroll, Phys. Rev. E 67, 066705 (2003); 
T. Ihle and D.M. Kroll, Phys. Rev. E 67, 066706 (2003). 
J. T. Padding and A. A. Louis, Phys. Rev. E. 74, 031402 
(2006). 

A. Malevanets and R. Kapral, J. Chem. Phys. 112, 7260 

(2000) . 

J. T. Padding and A. A. Louis, Phys. Rev. Lett. 93, 
220601 (2004); J.T. Padding and A. A. Louis, Phys. Rev. 
E. 77, 011402 (2008). 

M. Hecht, J. Harting, T. Ihle, and H. J. Herrmann, 
Phys. Rev. E 72, 011408 (2005). 

M. Ripoll, K. Mussawisade, R.G. Winkler, and G. Gomp- 

per, Phys. Rev. E 72, 016701 (2005). 

T.S. Squires and SR. Quake, Rev. Mod. Phys. 77, 977 

(2005) . 

E. S. Boek, J. Chin and P. V. Coveney, Int. Journ. Mod. 
Phys. B 17, 99 (2003). 

M. Venturoli, E. S. Boek, Elsevier Physica A 362 23, 

(2006) . 

G. A. Bird, Phys. Fluids 13, 2676 (1970). 



[24] G.A. Bird, Molecular Gas Dynamics and the Direct Sim- 
ulation of Gas Flows (Claredon Press, Oxford, 1994). 

[25] A. Succi, The Lattice Boltzmann Equation for Fluid Dy- 
namics and Beyond, Oxford University Press, Oxford 
(2001). 

[26] A.J.C. Ladd, Phys. Rev. Lett. 70, 1339 (1993) 
[27] A.J.C. Ladd and R. Verberg, J. Stat. Phys. 104, 1191 
(2001). 

[28] M.E. Cates et al, J. Phys.: Condens. Matter 16, S3903 
(2004); R. Adhikaril, K. Stratford, M. E. Cates, and A. 
J. Wagner, Europhys. Lett. 71, 473 (2005); K. Stratford 
et al, Science 309, 2198 (2005). 

[29] V. Lobaskin and B. Diinweg, New J. of Phys. 6, 54 

(2004) ; V. Lobaskin, B. Diinweg, and C. Holm, J. Phys.: 
Condens. Matt. 16, S4063 (2004). 

[30] F. Capuani, I. Pagonabarraga, and D. Frenkel, J. Chem. 

Phys. 121, 973 (2004). 
[31] A. Chatterji and J. Horbach, J. Chem. Phys. 122, 184903 

(2005) . 

[32] CP. Lowe, Europhys. Lett. 47, 145 (1999). 

[33] N. Kikuchi, C. M. Pooley, J. F. Ryder, and J. M. Yeo- 

mans, J. Chem. Phys. 119, 6388 (2003). 
[34] CM. Pooley and J.M. Yeomans, J. Phys. Chem. B 109, 

6505 (2005). 

[35] M. P. Allen and D. J. Tildesley, Computer simulation of 
liquids (Clarendon Press, Oxford, 1987) 

[36] J. T. Padding, A. Wysocki, H. Lowen, and A. A. Louis, 
J. Phys. -.Condens. Matter 17, S3393 (2005). 

[37] J. P. Hansen, I. R. McDonald Theory of Simple Liquids 
Academic Press Inc. (1986) 

[38] E. Tiizel, M. Strauss, T. Ihle, D.M. Kroll, Phys. Rev.E, 
68, 036701 (2003). 

[39] W. B. Russell, D. A. Saville, and W. R. Showalter, Col- 
loidal Dispersions (Cambridge University Press, Cam- 
bridge, England, 1989). 

[40] B.J. Alder and T.E. Wainwright, Phys. Rev. A 1, 18 
(1970). 

[41] M.H.J. Hagen, I. Pagonabarraga, CP. Lowe, D. Frenkel: 

Phys. Rev. Lett 78, 3785 (1997). 
[42] J. Sane, Phd Thesis (to be published) (2008). y 
[43] L.Bocquet, J-L. Barrat: J. Phys. Condens. Matter 8, 

9297 (1996) 

[44] P.M. Bungay and H. Brenner, Int. J. Multiph. Flow 1, 
25 (1973). 



